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Abstract 

Various authors have shown that, near the onset of a period-doubhng bifurcation, small perturbations 
in the control parameter may result in much larger disturbances in the response of the dynamical system. 
Such amplification of small signals can be measured by a gain defined as the magnitude of the disturbance 
' in the response divided by the perturbation amplitude. In this paper, the perturbed response is studied 

using normal forms based on the most general assumptions of iterated maps. Such an analysis provides 
a theoretical footing for previous experimental and numerical observations, such as the failure of linear 
analysis and the saturation of the gain. Qualitative as well as quantitative features of the gain are 
exhibited using selected models of cardiac dynamics. 

^ ! 1 Introduction 
O, 

^\ We study a dynamical system described by a one-parameter family of iterated maps 

o ■ 

^ : Xn+i = /(a;„;/i) , (1) 

where the state x may be one- or multi-dimensional. Suppose that a fixed point of the map experiences a 
supercritical period-doubling bifurcation when the control parameter fi passes through a critical value /ibif- 
^ ■ Near onset of the bifurcation, the system is sensitive to small alternating perturbations applied to /x [1-14]. 

^ I Under such alternating perturbations, the governing map becomes 

><■ x„+i=/(x„;/i+(-l)"5), (2) 

H ; 

where 6 is the perturbation amplitude. Consequently, the steady-state solution of the system alternates 
between two states, a^cvon and Xodd, in the following manner: 

a^ovon = / (a;odd; M ~ ^) (3) 
X odd = f {x even; ^J■ + S) . (4) 

The magnitude of deviation between Xeven and Xodd could be many times larger than the perturbation 
amplitude S, an effect which is known as pre-bifurcation amplification [7-9] . To characterize this amplification, 
one may define a gain using the ith component of x as follows 



*^cvcn 



dd 



26 



(5) 
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In this work, we explore qualitative changes in the gain F under changes in the control parameter and 
perturbation amplitude S. 

It has been shown that pre-bifurcation amplification is a universal phenomenon existing in discrete [1-6] 
as well as continuous [7-10] dynamical systems. Many authors utilized one-dimensional maps to illustrate 
the scaling of parameters. Such one-dimensional maps either take the most general form as in Heldstab et al. 
^ or take the normal form as in Bryant and Wiesenfeld 9 . Since the normal form in 9 was not related to 
the full-dimensional underlying model, it is not clear how the scaling of parameters is connected to physical 
properties of the system. 

Although the fundamental scaling law behind prebifurcation amplifiction has been shown in the literature, 
explicit solutions of the gain F are not known for arbitrarily-given iterated maps. In this paper, we aim to 
derive an analytical formula of the gain based on multi-dimensional iterated maps, requiring only minimum 
assumptions. Using numerical simulations as well as experimental work pi 1111 [T^ . a few researchers have 
shown that the correct scaling law must resort to a nonlinear analysis. We illustrate, for the first time, the 
failure of the linear analysis through an order analysis, which reveals that the range of validity for the linear 
analysis is extremely limited and dramatically shrinks to zero as the bifurcation point is approached. We 
also explain the saturation of amplification previously observed in |13[ [T^ using the developed scaling law. 

Another objective of the current work is to examine period-doubling bifurcations in cardiac dynamics, 
i.e., alternans in cardiac tissues, which are characterized by the long-short beat-to-beat alternation in the 
duration of cardiac action potentials [15-18]. Cardiac alternans have been recognized as a possible initiator of 
fatal arrhythmia, which is the number one cause of death in the United States 19' . We apply the theoretical 
results to two cardiac mappings to help gain an understanding of the origin and control of instabilities in 
cardiac tissues. 

The paper is organized as follows. In Section 2, we compute the response of alternating perturbations 
using linear analysis of the map. A nonlinear analysis based on higher-order approximation of the map is 
then performed in Section 3. A summary and discussion is presented in Section 4. 



2 Linear Analysis 



At fi — /ibif , denote the fixed point of the map by Xhn. Then the Jacobian fx [xhn] fJ-hn) has one eigenvalue 
equal to —1 with the corresponding right and left eigenvectors denoted by (t>i and ^i, respectively. Near 
onset of the period-doubling bifurcation, the response of the perturbed map |5J exhibits pre-bifurcation 
amplification. The following theorem computes the amplification gain based on linearization of the perturbed 
map. 



Theorem 1. In linear approximation, as fi 



^Mf, the gain defined in Eq. ^ diverges as 
k 



F : 



where the constant k is determined by data from the mapping. Specifically, 



k = 



where 



fxx {Xhif'i Mf)i/) ■ -qJ^ il^bif) + fxfi (Xbifi l^bif) : 



(6) 



(7) 



(8) 



(/i) is a fixed point of map 0). 



In Theorem 1, 0^*'' represents the ith component of 0i. For one-dimensional maps, the eigenvectors are 
01 = V'l = 1 and thus 

, _ fn (xhit; Mbif) 



fxx [Xhif, fJ-hif) ■ ^ (^bif) + /s/i {Xhif, Mbif) 



(9) 
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Note, here, subscripts indicate partial derivatives of the function /. This theorem indicates that the magni- 
tude of the gain depends on the properties of the system at the bifurcation point. Moreover, one can make 
two striking observations regarding the hnear approximation: i) the gain does not depend on the amphtude 
6 and ii) the gain approaches oo as fj, —>■ fibn- As shown at the end of this section, both observations compare 
poorly with numerical simulations. To rectify the discrepancy, we are led to the higher-order nonlinear 
analysis of the gain in Section 3. 



2.1 Proof of Theorem 1 

We prove the theorem by a perturbation-theory calculation. We assume the map ^ is A^-dimensional. 
Denote the fixed point of map as x, (/i), so that 

a;* (/") = / (2;* (a*) ; m) • (10) 

Let us assume the fixed point is stable when fi is greater than but close to /ibif (otherwise we could introduce 
a new control parameter a,s p, = —^i). We then study the response of the perturbed map |(5J. Under small 
perturbation, the response of the perturbed system is a small deviation (denoted by dn) from the unperturbed 
fixed point. Let us write a;„ a;, + which upon substitution into Eq. ((21, yields 

d„+i = L ■ d„ + (-1)" (x,;^) 5, (11) 

where 

L^f^{x.-li). (12) 

To account for the alternating feature of the perturbation, we decompose the steady state of dn into an 
oscillating component (—1)" v and a mean value m: i.e., dn — (—1)" v -\- m. Substituting dn into Eq. (jll|l 
yields 

m = L ■ m (13) 
V = -L-v- ff,{x;n)5. (14) 

It follows that m = in this linear approximation. 

In the following, we find a solution for v. Because fj, is close to /^bifj L can be approximated as 

L = Lo + (/i - /ibif) (15) 

where 

Lo = fx {xhii; A^bif) (16) 

Ll = fxx {Xhii; /^bif) • (Mbif) + fxt^ {Xhu; Mbif) • (17) 

Implicit differentiation of Eq. I|10(l yields 

(A^bif) = {I- .fx {xhiC, Mbif)) • (a;bif; Mbif) , (18) 

where I stands for the identity matrix. Denote the eigenvalues of by — 1 = Ai < A2 < ... < Ajy < 1 
and denote the corresponding right and left eigenvectors hy (pi, (j)2, (f>N and ipf, tp2, • ■ • , V'wj respectively. 
Assume the eigenvectors are normalized such that ipf ■ (f)j = 1 \i i = j and ipf ■ (pj = if i 7^ j. We can 
decompose v as 

N 

v^"^s,(l),. (19) 

1=1 

Substituting Eq. (HHJ into Eq. (|Tlj) yields 

N N 

^ (I + Lo) • (t>iSi + ip- Mbif) ^ Ll ■ (piSi = -ff^ {x; S. (20) 

i=l i=l 
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The dot product on both sides of Eq. H2U|I with tpj ioi j — 1, . . . , N yields 

TV 

(1 + Xj) Sj + (Ubif) ^ ipj ■ Li ■ (j)iSi = -V'J • {x; n) S. (21) 



Because 1 + Ai = 0, it follows that S2, ■ ■ •, sat are negligible compared to Si. To lowest order, solution of si 
is 

= i^I-Ui^;t^) (22) 

ilJ- - Mbif) {tPl ■ Li ■ (jji) 

Thus, an approximation of v is 

i'T ■ /m (a^; m) 



(/i - Mbif) ("01 ■ Li ■ (pi) 

In steady state, the solution of the perturbed system alternates between two states, Xcvon and Xodd, where 

■^cvcn X^ -\- 77, (^^) 

a^odd = a;* - w. (25) 



Recalling definition ^ of the gain yields 

r = 

where 







k 


s 




M - Mbif 


I ■ U (a^bif ; /^bif) 






■ 01 



(26) 
(27) 

2.2 An Example Illustrating the Limitation of the Linear Analysis 

As the first numerical example, we consider a one-dimensional cardiac map (16j : 

^„+i = A„ax-«e-^"/" (28) 

where £>„ = B — An- Here, An is the nth duration of action potential of a cardiac cell under an electrical 
pacing with period B (known as the basic cycle length); £)„ is the nth diastolic interval. Using the parameters 
Anax = 392.0 ms, a = 525.3 ms, and r = 40.0 ms, Hall and Gauthier showed that a period-doubling 
bifurcation occurs at B = B^if ^ 455 ms and the period-one fixed point of the system is stable for B > Bbif. 
We refer interested readers to ^5] for a bifurcation diagram exhibiting the period-doubling bifurcation and 
other properties of this map. 

It is our interest to study the response of this cardiac map when the basic cycle length B is perturbed 
with an alternating perturbation in such a way that D„ = B + (— 1)"5 — A„. Following Theorem 1 and 
evaluating the corresponding derivatives, one obtains the amplification gain as 

" (29) 



B - Bbif 

Again, the theorem predicts i) the gain does not depend on perturbation amplitude S and ii) the gain 
approaches oo when the control parameter B Bhu- For a given B, the magnitude of the gain depends on 
the value of the parameter r. To verify these observations, we compare results from theory and numerical 
simulations, as shown in Fig. ^ Here, Fig. ^ (a) shows variation of the gain under changes in the control 
parameter B when the perturbation amplitude S takes on various constant values and Fig. ^ (b) shows 
variation of the gain under changes in S when B takes on various constant values. In both figures, numerical 
simulations indicate that the gain strongly depends on S. It is also noted that the gain is not cx3 at i? = i?bif 
except when S — 0. These qualitative differences between the linear analysis and numerical simulations 
therefore demand a higher-order nonlinear analysis of the gain, which will be presented in the following 
section. 
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AB (ms) 5 (ms) 



Figure 1: Variation of the gain for the 1-D cardiac map as function of (a) AB — B — Bhn and (b) perturbation 
amphtude 5: the solid curves correspond to results from the linear analysis and the dots correspond to 
numerical simulations. Note that, in (b), the linear analysis predicts an infinitely large gain for Ai? = and 
any 5. 




AB (ms) 5 (ms) 

Figure 2: Variation of gain for the 1-D cardiac map as function of (a) Ai? = B — B\,i{ and (b) perturbation 
amplitude S: the solid curves correspond to results from the higher-order nonlinear analysis and the dots 
correspond to numerical simulations. 



3 Higher-Order Analysis 

Again, denote the fixed point of the map Q at = ^bif by Xhn- And again, denote by (f>i and if^i the right 
and left eigenvectors associated with the —1 eigenvalue of the Jacobian (a^bif ; Mbif ) • Near onset of the 
period-doubling bifurcation, the response of the perturbed map |(2J) exhibits pre-bifurcation amplification. 
The following theorem computes the amplification gain based on a higher-order nonlinear analysis of the 
perturbed map. 

Theorem 2. To leading order, the gain defined in Eq. |3|) satisfies the following relation, 

c6^r^ + {^i-^iMf}r-\k\ = o. (30) 

The coefficients in Eq. are defined in the following subsection. 

(i) 

Here, again, represents the ith component of (pi. For the special case of one-dimensional maps, the 
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eigenvectors are = ipi ^ 1 and thus 



fc = ^^ /1\ , . : (31) 



//i {Xhif, Mbif) 
fxx [Xhif, fJ-hif) ■ ^ (A'bif) + fxfi [xhif, Mbif) 

{l/2fxx {Xbn; fJ-hif)f + l/6fxxx {Xhif, Mbif) 



fxx (xhif, fJ-bif) • ifJ-bif) + /k^ [Xbif] fJ-hif) 



(32) 



Theorem 2 clearly indicates that the gain strongly depends on S as well as fi. Especially at the bifurcation 
point, where fi = fihif, the gain approaches to oo as a function 5~'^^^. As shown by numerical examples at 
the end of this section, the higher-order analysis agrees well with numerical simulations. 



3.1 Proof of Theorem 2 

Heldstab et al. ^ seemingly were the first to conduct a higher-order nonlinear analysis of the perturbed 
map 121) • Letting /i = /ibif, they calculated the response of alternating perturbations in one-dimensional 
maps. It was shown that deviation between Xovcn and Xodd is proportional to 5^^^. Following the same 
line of approach, we study the nonlinear response of maps of any dimensions and consider the influence of 
both the control parameter fj, and the perturbation amplitude S. Naturally, the first attempt of higher-order 
analysis is to include only quadratic terms, which, as shown in the following, modifies the expressions for 
the median of the response as well as modulates its amplitude. A retrospective examination reveals that the 
contribution to the response from the cubic term is comparable to that from the quadratic term and thus 
the cubic term has to be taken into account. 

Again, we assume the perturbed response is a deviation (i„ away from the unperturbed response x*, 
defined in Eq. H1U|I . Substituting Xn = x^ + d„ into Eq. ^ and keeping the expansion of / up to cubic 
order in dn yields 

dn+i ^ L ■ d„ + Q ■ dn ■ d„ + C ■ dn ■ dn ■ dn + (-1)" fp,{x; fi) S + . . . , (33) 

where 



L = fx {x; fi) 

Q = lfxx{x■,^l) (34) 

C = ^fxxx {x; n) . (35) 
6 

To account for the alternating feature of the perturbation, we decompose the steady state of dn into an 
oscillating component (—1)" v and a mean value m. Substituting dn = (—1)" v + m into Eq. (|33|l yields 

m = L- m-\-Q-v-v (36) 
V — —L ■ V — 2Q ■ m ■ V 

- C ■ V ■ V V - ff^{x;fj,)6. (37) 

Because ^ is close to /^bif, L can be approximated as 

i = io + (m - Mbif) -^1, (38) 

where Lq and Li are defined in Eqs. lfTB|l and l(T7|) . Again, denote the eigenvalues of Lq by — 1 = Ai < A2 < 

. . . < A AT < 1 and denote the corresponding right and left eigenvectors hy (j)i,(j>2, ■ ■ ■ ,<I>n and ipi, tp2, ■ ■ ■ , V'ati 

respectively. We assume the eigenvectors are normalized such that i/jf ■ (pj — 1 if i — j and ipf ■ — if 

It follows from Eq. (|36|l that, in steady state. 



Q-v-v, (39) 
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where 



Q = (I-Lo)"' -Q. (40) 



To solve for the steady-state solution of w, we again decompose v as 

N 



V = Ys^(|)^ = <^ ■ S, (41) 



i=l 



where $ — (01, . . . , ^jv) and s — {si,...sn) ■ Substituting Eq. H41|) into Eq. H37f) and taking the dot- 
product of this expression with tpj for j ^ 1, . . . , N on both sides yields 



N 



(1 + Aj) Sj + Aibif) V-J • Li ■ 0,s, (42) 
1=1 

= -2vj • g • (g • ($ • s) • ($ • s)) . ($ • s) (43) 

- • C • ($ . s) . ($ . s) . ($ . s) - • U (a^; ^- (44) 

Because 1 + Ai = 0, it foUows that §2, . . ., sat are higher order compared to si and thus v ~ 4>iSi. Then, 
the governing equation for si is, to lowest order, 

(g + c)sl + iti- fibit) Lsi + ^p'[-u f^)s = o, (45) 

where 

L ^ ■ Li ■ cj), (46) 

g = 2^-^ Q • (Q • 01 • '/'i) • ^1 (47) 
C = (48) 

The perturbed response governed by Eq. H45|) is valid for both supercritical and subcritical period- 
doubling bifurcations. However, for a supercritical period-doubling bifurcation such as the one mediating 
cardiac alternans, si = is the only solution when 6 = and /i > /ibif; that is to say, before the occurrence 
of a supercritical period-doubling bifurcation, the only solution of the unperturbed map is the period-one 
solution. Thus, for a supercritical period-doubling bifurcation, it follows from Eq. 1451) that 

{Q + C)L>0. (49) 

In the following, we restrict our attention to supercritical period-doubling bifurcations to simplify the result. 
In steady state, the solution of the perturbed system alternates between two states, Xeven and Xodd, where 

a;ovcn = a:* + w. (50) 
a^odd = - w. (51) 



Recalling the definition of the gain yields 



(52) 



Substituting Eq. into yields 

c52r3 + (^-Aibif)r- = 0, (53) 

where 



(54) 



and k is defined in Eq. (|27|l . 
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3.2 Examples Applying the Nonlinear Analysis 
3.2.1 One-Dimensional Cardiac Model 

Using the higher-order analysis, we reexamine the one-dimensional cardiac model described by Eq. H28|l . 
Results from the higher-order analysis agree well with numerical simulations, as seen in Fig. |21 It is observed 
from Eq. H53I) . as well as the numerical results, that the gain has a strong dependence on S and the gain 
in general is not oo at _B = By^n except when (5 = 0. The discrepancy between the linear and nonlinear 
analyses can be seen clearer in Fig. |31 which shows the error between the linear and nonlinear predictions in 
two-dimensional parameter space {AB,5). The figure shows that the error increases as 6 increases and for 
the same S the error decreases as AB increases. 




Figure 3: Discrepancy between the linear and nonlinear analyses for the the 1-D cardiac map. Here, the 
error is defined as (Fl — Fn) /Fn, where Fl is the result from the linear analysis and Fn is the result from 
the nonlinear analysis. 



3.2.2 Two-Dimensional Cardiac Model 

We now consider a two-dimensional map presented by Chialvo et al. IP and later studied by Hall and 
Gauthier 

An+i = (1 - M„+i) (A„ax - a e-^"/") , (55) 

M„+i = [1 - (1 - M„)e-4"/^^]e-^"/^^ (56) 

where Z)„ = B — An, An is the nth duration of action potential of a cardiac cell, B is the basic cycle 
length representing the period of applied electrical pacing, and Dn is the nth diastolic interval. Using the 
parameters A^ax — 490.9 ms, a — 569.0 ms, ti — 64.0 ms, T2 = 38.0 ms, I?min = 38.0 ms. Hall and Gauthier 
showed that a transition to alternans occurs at i? = i?bif ~ 455 ms ^Hl- We refer interested readers to 
[1,6' for a bifurcation diagram exhibiting the period-doubling bifurcation and other properties of this map. 
We study the response of this map under alternating perturbations. Figure 01 shows the variation of the 
amplification gain under changes in B and 6. Again, it is seen that good agreement is achieved between the 
higher-order analysis and the numerical simulations. 

4 Conclusion and Discussion 

Using both linear and higher-order analyses, we have investigated the pre-bifurcation amplification near a 
period-doubling bifurcation in multi-dimensional iterated maps. Through numerical examples, it has been 
estabilished that the qualitative features of the amplification gain can only be captured using the higher- 
order analysis. In general, the gain is a nonlinear function of fi and S as seen from Eq. (|53f) . In the {fJ-,5) 
parameter region, the gain is finite everywhere except at (/ibif: 0)- Moreover, the rate of divergence as (/i, S) 
approaches to this singular point depends on the path taken. Under two situations, this nonlinear function 
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Figure 4: Variation of the gain for the two-dimensional cardiac map as function of (a) AB = B — Bbif and 
(b) perturbation amphtude S: the sohd curves correspond to results from the higher-order nonlinear analysis 
and the dots correspond to numerical simulations. 



can be significantly simplified. First, when |/i — /ibif| '3> S^^^, the first term in Eq. (|53|l becomes negligible 
compared to the others, resulting in the simplified expression 

(A*-/ibif)r-|fc| =0. (57) 

It follows that the gain can be approximated by 

r = ^— , (58) 

fJ. ~ Mbif 

which is the result of linear analysis obtained in Section 2. On the other hand, when |/i — /ibifl ^ S^^^, the 



Figure 5: Range of validity: the linear analysis is valid in the hatched area; the higher-order analysis is valid 
in the whole region. Here, it has been assumed that region of interest is when fi > /ibif. 

second term in Eq. (|53|l becomes negligible compared to the others, resulting in the simplified expression 

c5^T^-\k\^Q. (59) 
It follows that the gain can be approximated by 




A special case of the latter situation is when /i = /ibif, which was studied by Heldstab et al. using one- 
dimensional maps T. Note that Heldstab et al. did not include the contribution due to the cubic order 
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expansion of v {i.e. the contribution of C), although their result is still valid for the logistic map because 
C = for this map. Our analysis allows us to predict the range of validity for the linear and the higher-order 
analyses, as shown in Fig. |S1 Especially, the linear analysis is only valid in a restricted area, which shrinks 
to zero as ^ approaches /^bif, cf. the numerical illustration in Fig. |21 

The two aforementioned special situations and the analysis of the range of validity lead to an interesting 
observation on the relationship between the gain F and the parameters. First, when 5 takes on a constant 
value, one can see from Figure El (a) that logF is almost proportional to log (/i — /^bif ) ^ for large ^i, and 
logF increases and saturates to a constant as /i decreases to a sufficiently small value. Conversely, when 
/i takes on a constant value, one can see from Figure El (b) that logF is almost proportional to log (5^^/"^ 
for large (5, and logF increases and saturates to a constant as 5 decreases to a sufficiently small value. We 
note that a similar saturation effect in pre-bifurcation amplification was previously observed by Kravtsov 
and Surovyatkina 13 and Surovyatkina .14. . who considered random noise in a one-dimensional map as /i 
is varied to approach /Zbif. Although the saturation of the gain is a universal property of dynamical systems, 
we note that it is very difficult to observe experimentally such phenomenon on paced cardiac tissues because 
of the existence of noise and the limitation on measurement resolution. 
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Figure 6: Saturation of the amplification gain in the two-dimensional cardiac map illustrated by (a) logF 
vs. log Ai3, where Ai? = B — Bhu and (b) logF vs. log a: the solid curves correspond to results from the 
higher-order analysis and the dots correspond to numerical simulations. 

Because pre-bifurcation amplification is a property of the period-doubling bifurcation, its observable 
region in terms of the control parameter and the perturbation amplitude is problem-dependent and limited 
in parameter space to a vicinity of the bifurcation point. Under some circumstances, de-amplification gains 
(gains less than 1) may be obtained for a control parameter /i sufficiently far away from /ibif. Moreover, 
other responses rather than the presumed period-two solutions may also exist. For example, as shown in 
Fig. [7| a period-doubling cascade to chaos is observed in the two-dimensional cardiac model under increase 
of the perturbation amplitude S. These observations suggest that caution should be taken when using 
pre-bifurcation amplification as a technique to detect the existence of period-doubling bifurcations. 
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